PROBLEMA - SERIE GAS
library(car)
library(urca)
library(forecast)
library(tseries)
library(ggfortify)
library(TSstudio)
library(highcharter)Producción de gas mensual en Australia 1956-1995
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1646 2675 16788 21415 38629 66600
Transformación de los datos
## [1] 0.08262296
## [1] 476
Identificación
par(mfrow=c(3,3))
plot(z1,type="o")
acf(z1, lag.max=40)
pacf(z1, lag.max=40)
plot(diff(z1),type="o")
abline(h=2*sqrt(var(diff(z1))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(z1))),col="red",lty=2)
acf(diff(z1), lag.max=40)
pacf(diff(z1), lag.max=40)
plot(diff(z1,12),type="o")
acf(diff(z1,12), lag.max=40)
pacf(diff(z1,12), lag.max=40)Test DF sobre la serie desestacionalizada
Dickey-Fuller
##
## Augmented Dickey-Fuller Test
##
## data: z1
## Dickey-Fuller = -0.72371, Lag order = 7, p-value = 0.9683
## alternative hypothesis: stationary
## Warning in pp.test(z1): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: z1
## Dickey-Fuller Z(alpha) = -29.784, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary
- De acuerdo a la prueba de Dickey- Fuller el valor p obtenido fue 0.9586 por lo que no es posible rechazar la hipotesis nula y diremos no es estacionaria
Dickey Fuller con una diferencia
## Warning in adf.test(diff(z1)): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff(z1)
## Dickey-Fuller = -17.802, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
## Warning in pp.test(diff(z1)): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: diff(z1)
## Dickey-Fuller Z(alpha) = -302.63, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary
Los resultados de la serie anterior nos llevan a concluir que las serie con una diferencia ordinaria es estacionaria
D-F y PP sobre la serie desestacionalizada
##
## Augmented Dickey-Fuller Test
##
## data: diff(z1, 12)
## Dickey-Fuller = -3.5139, Lag order = 7, p-value = 0.04111
## alternative hypothesis: stationary
## Warning in pp.test(diff(z1, 12)): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: diff(z1, 12)
## Dickey-Fuller Z(alpha) = -51.992, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary
Las pruebas de Dickey Fuller y Phillips Perron nos indican que la serie z1 con una diferencia estacional es estacionaria
## Warning in adf.test(diff(diff(z1, 12))): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff(diff(z1, 12))
## Dickey-Fuller = -7.8325, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
par(mfrow=c(4,3))
plot(z1,type="o")
acf(z1, lag.max=40)
pacf(z1, lag.max=40)
plot(diff(z1),type="o")
abline(h=2*sqrt(var(diff(z1))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(z1))),col="red",lty=2)
acf(diff(z1), lag.max=40)
pacf(diff(z1), lag.max=40)
plot(diff(z1,12),type="o")
acf(diff(z1,12), lag.max=40)
pacf(diff(z1,12), lag.max=40)
plot(diff(diff(z1,12)),type="o")
abline(h=2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
acf(diff(diff(z1,12)), lag.max=40)
pacf(diff(diff(z1,12)), lag.max=40)par(mfrow=c(1,3))
plot(diff(diff(z1,12)),type="o")
abline(h=2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(diff(z1,12)))),col="red",lty=2)
acf(diff(diff(z1,12)), lag.max=40)
pacf(diff(diff(z1,12)), lag.max=40)## Series: z1
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.3752 -0.8589
## s.e. 0.0450 0.0451
##
## sigma^2 = 0.0122: log likelihood = 356.1
## AIC=-706.2 AICc=-706.14 BIC=-693.78
SARIMA(0,1,1)(0,1,1) BIC -347.7133 Lags:1,1
SARIMA(0,1,5)(0,1,1) BIC -329.9128 Lags:1:5,1
SARIMA(0,1,5)(0,1,1) BIC -341.6644 Lags:1,3,5;1
SARIMA(0,1,5)(0,1,1) BIC -345.9592 Lags:1,5;1
SARIMA(0,1,11)(0,1,1) BIC -344.8457 Lags:1,5, 10, 11;1
SARIMA(0,1,11)(0,1,1) BIC -349.3178 Lags:1,5, 11;1
SARIMA(0,1,16)(0,1,1) BIC -344.7575 Lags:1,5,11,16;1
SARIMA(16,1,11)(0,1,1) BIC -347.3572 Lags:16;1,5,11;1
SARIMA(16,1,19)(0,1,1) BIC -342.4163 Lags:16;1,5,11,18,19;1
SARIMA(19,1,18)(0,1,1) BIC -342.1835 Lags:16,19;1,5,11,18;1
SARIMA(16,1,18)(0,1,1) BIC -347.9725 Lags:16;1,5,11,18;1
Ajuste del Modelo
#modelo1<-stats::arima(z1,
#order=c(16,1,18),
#seasonal=list(order=c(0,1,1),
# period=12),
#fixed=c(rep(0,15),NA,NA,0,0,0,NA,0,0,0,0,0,NA,0,0,0,0,0,0,NA,NA))
#modelo1
#tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
#1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))
#BIC(modelo1)modelo1<-stats::arima(z1,
order=c(0,1,11),
seasonal=list(order=c(0,1,1),
period=12),
fixed=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA))
modelo1##
## Call:
## stats::arima(x = z1, order = c(0, 1, 11), seasonal = list(order = c(0, 1, 1),
## period = 12), fixed = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8 ma9
## -0.3891 0.0127 -0.0430 0.0371 0.0964 -0.0554 0.0552 0.0282 0.0585
## s.e. 0.0460 0.0519 0.0511 0.0519 0.0510 0.0518 0.0550 0.0524 0.0615
## ma10 ma11 sma1
## 0.0588 0.1063 -0.8965
## s.e. 0.0540 0.0592 0.0392
##
## sigma^2 estimated as 0.01149: log likelihood = 367.3, aic = -708.59
tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))## ma1 ma2 ma3 ma4 ma5 ma6
## 2.220446e-16 4.033875e-01 2.006337e-01 2.373491e-01 2.977826e-02 1.427718e-01
## ma7 ma8 ma9 ma10 ma11 sma1
## 1.580442e-01 2.954415e-01 1.709673e-01 1.383101e-01 3.667343e-02 0.000000e+00
## [1] -654.8039
Test de Autocorrelacion de Ljung-Box
\(H_0\): \(r_1=r_2=r_3=...=r_{lag}=0\)
\(H_a\): Al menos una es diferente de cero
#autoplot(modelo1)
#LBQPlot(et)
##
## Box-Ljung test
##
## data: et
## X-squared = 0.92628, df = 11, p-value = 1
Test de Normalidad basado en Sesgo y Curtosis
\(H_0=\): Datos provienen de una Dist. Normal
\(H_a\): Los datos no provienen de una Dist. Normal
##
## Jarque Bera Test
##
## data: et
## X-squared = 595.33, df = 2, p-value < 2.2e-16
Test de Aleatoriedad
\(H_0:\) Residuales exhiben un comport. de Aleatoriedad \(H_a:\) Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)
##
## Runs Test
##
## data: as.factor(sign(et))
## Standard Normal = 0.64872, p-value = 0.5165
## alternative hypothesis: two.sided
inversa de boxcox
lambda <- 0.082
predt <- forecast(modelo1,h=12)
z2inv <- forecast::InvBoxCox(predt$mean,lambda)
z2inv.li <- forecast::InvBoxCox(predt$lower[,2],lambda) #Linf IC 95%
z2inv.ls <- forecast::InvBoxCox(predt$upper[,2],lambda) #lsUP ic 95%
z1.fit <- forecast::InvBoxCox(x1.fit,lambda)
plot(ts(c(gas,z2inv),start=c(1956,1),freq=12), type="l", col="blue", lwd=2,
main="Pronóstico h=12 Pasos al Frente Gas",
xlab="Anual",
ylab="",
ylim=c(min(gas,z2inv.li,z2inv.ls),max(gas,z2inv.li,z2inv.ls)))
polygon(c(time(z2inv.li),rev(time(z2inv.ls))),
c(z2inv.li,rev(z2inv.ls)),
col="gray", border=NA)
lines(z2inv, type="b", col="blue", lwd=2)
lines(z1.fit, type="l", col="red", lty=2, lwd=3) ## Jan Feb Mar Apr May Jun Jul Aug
## 1995
## 1996 41845.87 43321.40 46885.40 48679.72 57202.42 61757.94 65369.27 63345.99
## Sep Oct Nov Dec
## 1995 55161.78 51423.89 47349.65 43455.69
## 1996
## Jan Feb Mar Apr May Jun Jul Aug
## 1995
## 1996 48014.92 50356.73 55016.67 57747.27 68477.69 74758.82 80105.46 78877.02
## Sep Oct Nov Dec
## 1995 60087.85 56874.84 53111.61 49299.27
## 1996
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## m s df
## -3.300286e-05 7.245611e-02 3.533389e+00
## ( 3.976008e-03) ( 4.464714e-03) ( 6.517475e-01)
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in dt((x - m)/s, df, log = TRUE): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## Warning in log(s): Se han producido NaNs
## m s df
## -3.300286e-05 7.245611e-02 3.533389e+00
## ( 3.976008e-03) ( 4.464714e-03) ( 6.517475e-01)
## [1] 447 170